////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//  ----------------------------- Ebrahim Foulaadvand, 29 July 2012 ------------------------------------- //
//                                                                                                        //
// The programme "DrivenDampedOscillator" evaluates the temporal evolution of a sinusoidally driven       //
// damped harmonic oscillator by two algorithms Euler-Richardson (ER) and 2nd order Runge-Kutta (RK2).    //
// It also computes the oscillator amplitude dependence on the driving force frequency.                   //
//                                                                                                        //                                                                                               //
//                                                                                                        //                                                                                                         //
////////////////////////////////////////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>

using namespace std;


double a (double &t,double &x,double &v); // a is the accelration function

double  tau=0.01,T=400,ac,k=9.,m=1.,w0=sqrt(k/m),x0=0.,v0=0.,gamma=0.5,beta=gamma/(2*m),omega,delta,
        A0=1.,E,kx1,kx2,lx1,lx2,tmid,xmid,vmid,Aanalytic,Anumeric=0.,omegamax=4.,delomega=0.1;

// x0=initial position, v0=initial velocity, tau=timestep, T=simulation time, E=total energy,
// m=oscillator mass, k=spring constant, gamma=damping coefficient, omega=frequency of external force.
// A0=F0/m where F0 is the driving force amplitude, Aanalytic=analytical amplitude, Anumeric=computed
// amplitude, omegamax=the maximum frequncy for which the amplitude is comuputed, delomega=omega increment.

int N=(T/tau), Nomega=int(omegamax/delomega); //N=number of time steps

main()

{

vector<double> v(N+1,0),x(N+1,0),Acomp(Nomega,0),Aanal(Nomega,0);

// arrays x and v store oscillator position and velocity. Array Acomp stores the amplitude for a given omega.

ofstream xER("xER.plt");    // output file of oscillator poistion (ER algorithm).
ofstream vER("vER.plt");    // output file of oscillator velocity (ER algorithm).
ofstream xRK2("xRK2.plt");  // output file of oscillator poistion (RK2 algorithm).
ofstream vRK2("vRK2.plt");  // output file of oscillator velocity (RK2 algorithm).
ofstream Energy("Energy.plt");
ofstream phase("phase-space.plt");
ofstream Ampcomp("Acomp-w gamma=0.5.plt"); // output file of computed oscillator amplitude vs omega.
ofstream Ampanal("Aanal-w gamma=0.5.plt"); // output file of analytic oscillator amplitude vs omega.


for(int s=0;s<Nomega;s++){ cout<<"s= "<<s<<endl;//getch();

omega+=delomega;

Aanalytic=A0/(sqrt((w0*w0-omega*omega)*(w0*w0-omega*omega) + 4*beta*beta*omega*omega));

//cout<<"Aanalytic= "<<Aanalytic<<endl;

Ampanal<<omega<<" "<<Aanalytic<<endl;


//------------------------ Euler-Richardson Method------------------------

x[0]=x0; v[0]=v0; Anumeric=0;

for(int t=0;t<N;t++){

tmid=t*tau + 0.5*tau;
xmid=x[t] + 0.5*v[t]*tau;
vmid=v[t] + 0.5*tau*a(t*tau,x[t],v[t]);

v[t+1]=v[t] + tau*a(tmid,xmid,vmid);
x[t+1]=x[t] + tau*vmid;


xER<<t*tau<<" "<<x[t]<<endl;
vER<<t*tau<<" "<<v[t]<<endl;
phase<<x[t]<<"  "<<v[t]<<endl;

if (t*tau>(10/beta) && t*tau<(10/beta) + 2*2*(M_PI/omega) && x[t]>0 && x[t+1]>x[t])  Anumeric=x[t+1];


}//end t

Ampcomp<<omega<<" "<<Anumeric<<endl; //cout<<"Anumeric ER= "<<Anumeric<<endl;


//------------------------- Runge-Kutta 2nd order Method ------------------------

x[0]=x0; v[0]=v0;

for(int t=0;t<N;t++){

kx1=tau*v[t];
lx1=tau*a(t*tau,x[t],v[t]);

kx2=tau*(v[t]+lx1);
lx2=tau*a(t*tau+tau,x[t]+kx1,v[t]+lx1);

x[t+1]=x[t]+0.5*(kx1+kx2);
v[t+1]=v[t]+0.5*(lx1+lx2);

xRK2<<t*tau<<" "<<x[t]<<endl;
vRK2<<t*tau<<" "<<v[t]<<endl;

if (t*tau>(10/beta) && t*tau<(10/beta) + 2*(M_PI/omega) && x[t]>0 && x[t+1]>x[t])  Anumeric=x[t+1];

}//end t

//Ampcomp<<omega<<" "<<Anumeric<<endl;

//cout<<"Anumeric RK2= "<<Anumeric<<endl;getch();

}// end s

cout<<"end."<<endl;getch();

return 0;

}

//----------------------------------------------------------------

double a (double &t,double &x,double &v) {

return ac=-k*x/m-(gamma/m)*v + A0*cos(omega*t);

}